cd C:\Users\crystal.zhan\Dropbox\Sheldon-Zhan\Homeownership\Data
* Use PUMA if county is unidentifiable
* Use the new aggregate disaster data that exclude disaster types that may not damage properties
* Remove heatwave from google trend search
* Cluster at the origin state level

set matsize 10000
capture log close
log using ownshp_disater0618_googletrend.log, replace
* Using the updated disaster data
* drop climnate zone in nnmatch

***Prepare Data***

use "C:\Users\crystal.zhan\Dropbox\Sheldon-Zhan\Migration\Data\acs0514_13.dta", clear
keep if migrate1d == 32

sort year statefip countyfips puma

merge n:1 year statefip countyfips puma using acs_hs0514_13.dta
drop _merge

* Identified vs unidentified counties
gen ifcnty = (countyfips ~= 0)
gen puma_cnty = puma if ifcnty == 0
replace puma_cnty = countyfips if ifcnty == 1

sort statefip ifcnty puma_cnty year
merge n:1 statefip ifcnty puma_cnty year using Disaster_puma_2yrlag0613.dta
drop if _merge == 2

local dsvar "damage_property damage_crop damage_total INJURIES_DIRECT INJURIES_INDIRECT DEATHS_DIRECT DEATHS_INDIRECT n_disa"
foreach var of varlist `dsvar'{
replace `var' = 0 if _merge == 1
}
drop _merge


* Labor market
merge n:1 statefip year using acs_state_labor_market.dta
drop if _merge == 2
drop _merge

* County latitude and longitude
sort statefip countyfips
merge n:1 statefip countyfips using county_lat_lon_13.dta
drop if _merge == 2
drop _merge
* Only the unidentified county cannot be merged

* Puma latitude and longitude
sort statefip puma
merge n:1 statefip puma using puma_lat_lon_13.dta
drop if _merge == 2
drop _merge

* Use the latitude & longitude or elevation of puma or county whichever is smaller
egen t = sd(puma_lat) if cnty_lat ~=. & puma_lat ~= ., by(statefip countyfips year)
gen latitude = cnty_lat
gen longitude = cnty_lon
replace latitude = puma_lat if cnty_lat == . | (t ~= . & t ~= 0)
replace longitude = puma_lon if cnty_lat == . | (t ~= . & t ~= 0)
drop t

*******
* County elevation
sort statefip countyfips
merge n:1 statefip countyfips using Avg_elevation_county_13.dta
drop if _merge == 2
drop _merge
* Only the unidentified county cannot be merged

* Puma elevation
sort statefip puma year
merge n:1 statefip puma year using Avg_elevation_puma_13.dta
drop if _merge == 2
drop _merge

* Use the elevation of puma or county whichever is smaller
egen t = sd(elev_in_ft_puma) if elev_in_ft ~=. & elev_in_ft_puma ~= ., by(statefip countyfips year)
gen elevation = elev_in_ft
replace elevation = elev_in_ft_puma if elev_in_ft == . | (t ~= . & t ~= 0)

drop t

***********

* Climate zone
sort statefip countyfips
merge n:1 statefip countyfips using climatezone0826_13.dta
drop if _merge == 2
drop _merge

sort statefip puma year
merge n:1 statefip puma year using climatezone_puma_13.dta
drop if _merge == 2
drop _merge

replace baclimatezone1 = baclimatezone_puma if baclimatezone1 == .
replace ieccmoistureregime1 = ieccmoistureregime_puma if ieccmoistureregime1 == .
replace ieccclimatezone = ieccclimatezone_puma if ieccclimatezone == .

drop baclimatezone_puma ieccclimatezone_puma ieccmoistureregime_puma baclimatezone ieccmoistureregime

*****************
* DMA
*****************
merge n:1 statefip puma using PUMA10toDMA95.dta
drop if _merge == 2
drop _merge

merge n:1 statefip countyfips using countytoDMA0226.dta
drop if _merge == 2
drop _merge

replace dma = dmaindex if dma == .

******
* Combine metropolitan areas and DMA
gen t = met2013
replace t = -dma if met2013 == 0 & dma ~= .
replace t = -statefip * 1000 if met2013 == 0
egen met_dma = group(t)
drop t

* Metropolitan & state
gen t = met2013
replace t = -statefip if met2013 == 0
egen met_stt = group(t)

*********

gen ownshp = (ownershp == 1)
replace ownshp = . if ownershp == .

egen stt_cnty = group(statefip ifcnty puma_cnty)
drop if stt_cnty == .

***Define Treatment

gen fatalities = DEATHS_DIRECT + DEATHS_INDIRECT

gen treat = 0
tab fatalities
replace treat = 1 if fatalities > 3

*************
* Prepare the control variables
tab hhtype, gen(hhtp)
drop hhtp1

* Education
gen yschl = 0 if educd <= 12
replace yschl = 2 if educd == 10
replace yschl = 1 if educd == 14
replace yschl = 2 if educd == 15
replace yschl = 3 if educd == 16
replace yschl = 4 if educd == 17
replace yschl = 5.5 if educd == 21
replace yschl = 5 if educd == 22
replace yschl = 6 if educd == 23
replace yschl = 7.5 if educd == 24
replace yschl = 7 if educd == 25
replace yschl = 8 if educd == 26
replace yschl = 9 if educd == 30
replace yschl = 10 if educd == 40
replace yschl = 11 if educd == 50
replace yschl = 12 if educd >= 60 & educd <= 64
replace yschl = 12.5 if educd == 65
replace yschl = 13 if educd == 71
replace yschl = 14 if educd == 81 | educd == 82 | educd == 83
replace yschl = 16 if educd == 101
replace yschl = 18 if educd == 114
replace yschl = 20 if educd == 115
replace yschl = 21 if educd == 116
label var yschl "Years of Schooling"

* Race
gen race_white = 0
replace race_white = 1 if race == 1 & hispan == 0
gen race_black = 0
replace race_black = 1 if race == 2 & hispan == 0
gen race_asian = 0
replace race_asian = 1 if (race == 4 | race == 5 | race == 6) & hispan == 0
gen race_other = 0
replace race_other = 1 if (race == 3 | race > 6) & hispan == 0
gen hisp = (hispan > 0)

tab baclimatezone1, gen(czone)
drop czone1

gen ifchld18 = (n_schlage > 0 | nchlt5 > 0)
label var ifchld18 "If there are children under 18"

tab metro, gen(mts)
drop mts1 mts2

* Generalize the metropolitan status
gen metro1 = (metro > 1 & metro ~= .)
replace metro1 = . if metro == .

gen immig = (citizen == 2 | citizen == 3)

**********************************
* Adjust for inflation
*******
gen cpi = 195.3 if year ==	2005
replace cpi = 201.6 if year == 2006
replace cpi = 207.3 if year == 2007
replace cpi = 215.303 if year == 2008
replace cpi = 214.537 if year == 2009
replace cpi = 218.056 if year == 2010
replace cpi = 224.939 if year == 2011
replace cpi = 229.594 if year == 2012
replace cpi = 232.957 if year == 2013
replace cpi = 236.736 if year == 2014

replace propinsr = propinsr/cpi * 195.3
replace rentgrs = rentgrs/cpi * 195.3
replace valueh = valueh/cpi * 195.3
replace owncost = owncost/cpi * 195.3
replace hhincome = hhincome/cpi * 195.3

* Property tax
gen l_m_proptx1 = round(l_m_proptx99_o)
gen l_m_proptx2 = 0 if l_m_proptx1 == 1
replace l_m_proptx2 = (l_m_proptx1 - 1)*50 - 25 if l_m_proptx1 < 22 & l_m_proptx1 > 1
replace l_m_proptx2 = (l_m_proptx1 - 22)*100 + 1050 if l_m_proptx1 >= 22 & l_m_proptx1 < 57
replace l_m_proptx2 = (l_m_proptx1 - 23)*100 + 1050 if l_m_proptx1 > 57 & l_m_proptx1 < 63
replace l_m_proptx2 = 5250 if l_m_proptx1 == 63
replace l_m_proptx2 = 5750 if l_m_proptx1 == 64
replace l_m_proptx2 = (l_m_proptx1 - 65)*1000 + 6500 if l_m_proptx1 > 64
replace l_m_proptx2 = 10000 if l_m_proptx1 == 69

foreach x in l_m_propinsr_o l_m_rentgrs_r l_m_proptx2{
replace `x' = `x'/1000
}


*Merge in Google Trends Data
cd "GoogleTrends"

gen disaster_year = year - 2 

sort migplac1 disaster_year
merge m:1 migplac1 disaster_year using google_allDisasters.dta
drop if _merge == 2

gen info_high1 = (blizzard>95 | flashflood>95 | flood>95 | highwinds>95 | hurricane>95 | icestorm>95 | tornado>95 | tropicalstorm>95 | wildfire>95 | winterstorm>95) 
gen info_high2 = (blizzard>90 | flashflood>90 | flood>90 | highwinds>90 | hurricane>90 | icestorm>90 | tornado>90 | tropicalstorm>90 | wildfire>90 | winterstorm>90) 
gen info_high3 = (blizzard>80 | flashflood>80 | flood>80 | highwinds>80 | hurricane>80 | icestorm>80 | tornado>80 | tropicalstorm>80 | wildfire>80 | winterstorm>80)

forvalues i = 1/3{
replace info_high`i' = . if _merge ~= 3
gen treat_info`i' = treat * info_high`i'
}


cd C:\Users\crystal.zhan\Dropbox\Sheldon-Zhan\Homeownership\Data\Results\20180618main

qui reg ownshp treat age yschl hhtp* hisp race_other race_asian race_black hhincome famsiz n_worker ifchld18 immig mts* l_m_propinsr_o l_m_rentgrs_r l_m_proptx2 l_p_farm i.year if year < 2013 & gq <= 2 & migrate1d == 32 & stt_cnty ~= . [aweight = hhwt]
gen insample0 = e(sample)


forvalues i = 1/3{

areg ownshp treat treat_info`i' age yschl hhtp* hisp race_other race_asian race_black hhincome famsiz n_worker ifchld18 immig mts* l_m_propinsr_o l_m_rentgrs_r l_m_proptx2 l_p_farm i.year if year < 2013 & gq <= 2 & migrate1d == 32 [aweight = hhwt], cluster(migplac1) absorb(stt_cnty)
qui outreg2 using google, keep(treat treat_info`i') tex(frag pretty) bdec(3) bfmt(f)  sdec(3) sfmt(f) symbol(***, **, *) bracket se addstat(Adjusted R-squared, e(r2_a)) 
test treat + treat_info`i' == 0
}



** Nearest neighbor matching
****
teffects nnmatch (ownshp year latitude longitude elevation age yschl hhincome ifchld18 famsize l_m_rentgrs_r l_m_proptx2 l_pop_growth l_p_unemp metro1) (treat) if year < 2013 & gq <= 2 & migrate1d == 32 [fweight = hhwt], vce(robust) ematch(year metro1) nneighbor(1) generate(match) osample(nomatch)
predict nn_dist*, distance

rename match1 omatch1
rename nn_dist1 onn_dist1

drop match* nn_dist*

egen t1 = rank( onn_dist1 ) if onn_dist1 ~= .
egen t2 = count( onn_dist1 ) if onn_dist1 ~= .
* Percentile of distance
gen pc_dist = t1/t2

**
* Trim sample using the same cutoffs as the baseline sample

forvalues i = 1/3{
display `i'
areg ownshp treat treat_info`i' age yschl hhtp* hisp race_other race_asian race_black hhincome famsiz n_worker ifchld18 immig mts* l_m_propinsr_o l_m_rentgrs_r l_m_proptx2 l_p_farm i.year if insample0 == 1 & onn_dist1 < 2.081718 [aweight = hhwt], cluster(migplac1) absorb(stt_cnty)
qui outreg2 using google, tex(frag pretty) keep(treat treat_info`i') bdec(3) bfmt(f) sdec(3) sfmt(f) symbol(***, **, *) bracket se addstat(Adjusted R-squared, e(r2_a)) append
test treat + treat_info`i' == 0

areg ownshp treat treat_info`i' age yschl hhtp* hisp race_other race_asian race_black hhincome famsiz n_worker ifchld18 immig mts* l_m_propinsr_o l_m_rentgrs_r l_m_proptx2 l_p_farm i.year if insample0 == 1 & onn_dist1 < 1.505425 [aweight = hhwt], cluster(migplac1) absorb(stt_cnty)
qui outreg2 using google, tex(frag pretty) keep(treat treat_info`i') bdec(3) bfmt(f) sdec(3) sfmt(f) symbol(***, **, *) bracket se addstat(Adjusted R-squared, e(r2_a)) append
test treat + treat_info`i' == 0

areg ownshp treat treat_info`i' age yschl hhtp* hisp race_other race_asian race_black hhincome famsiz n_worker ifchld18 immig mts* l_m_propinsr_o l_m_rentgrs_r l_m_proptx2 l_p_farm i.year if insample0 == 1 & onn_dist1 < 1.217786 [aweight = hhwt], cluster(migplac1) absorb(stt_cnty)
qui outreg2 using google, tex(frag pretty) keep(treat treat_info`i') bdec(3) bfmt(f) sdec(3) sfmt(f) symbol(***, **, *) bracket se addstat(Adjusted R-squared, e(r2_a)) append
test treat + treat_info`i' == 0
}

log close
